clear;

psi=1.5;
gamma=8;
rho=0.042;
sigma=0.09;
etaI=5;
etaX=5;

beta0=39;
deltaI=0.06;
deltaX=0.06;

A=0.26;
%lal0=1/29.5;
beta1=0;
%=(1+beta0)*(1-lal0)/(A*lal0);

ds=0.00005;
ss=ds:ds:0.1;
cbar=0.04;
lambda0G=0.05;
lambda0B=2;



lambda0=lambda0B;
zeta0=1;
zeta1=0.3;





for j=1:length(ss)
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
lambda(j)=lambdas;
rp=gamma*sigma^2+lambdas*(betas/(betas-gamma)-betas/(betas+1-gamma))-lambdas/(betas+1);
computestarB;
ixx(j)=ix;
qx(j)=1/(1-etaI*ix);
xx(j)=x;
bx(j)=(A-ix-x)^(1/(1-psi))*(rho*(1/(1-etaI*ix)))^(-psi/(1-psi));
dbx(j)=(A-ix-x)^(-1/psi)*rho/(1-etaX*x/s)*bx(j)^(1/psi);





phii=ix-etaI*ix^2/2-deltaI;
 gx(j)=phii-lambdas/(betas+1);
 cx(j)=(bx(j)^(1-1/psi)/rho*(1-etaI*ix))^(-psi);

 
    r=cx(j)/qx(j)+gx(j)-rp;
  theta=x*(1-etaI*ix);
  qB(j)=(A-ix)/(r+rp+theta-phii+lambdas/(betas+1));
  qG(j)=(A-ix-x)/(r+rp-phii+lambdas/(betas+1));
  qs(j)=1/(1-etaI*ix);
 rpx(j)=rp;
 rpB(j)=rp+theta;
rm(j)=r+rp;
 ux(j)=((r+rp-phii+lambdas/(betas+1))*rho^(-psi))^(1/(1-psi));
 %cx(j)=(r+rp-phii+lambda/(betas+1))*qB(j);
 alphastarx(j)=x/(A-ix-cbar);


end
bs=bx;
dbs=dbx;

bs0=bx;
ixx0=ixx;
xx0=xx;

vol=sqrt(sigma^2+2*lambdas/betas^2)


for j=1:length(ss)-1
dbxnu(j)=(bx(j+1)-bx(j))/ds;
end
dbxnu(j+1)=2*dbxnu(j)-dbxnu(j-1);
i=1;
while dbxnu(i)>dbx(i) && i<length(ss)
i=i+1;
end
flagstar=i;

xx(flagstar)
ss(flagstar)

for j=1:length(ss)
 phiI(j)=ixx(j)-etaI/2*ixx(j)^2-deltaI;
 phiX(j)=xx(j)/ss(j)-etaX/2*xx(j)^2/ss(j)^2-deltaX;
end






subplot(3,2,1)
plot(ss,ixx,'LineWidth',1.5)
hold on;
xlabel('s')
title('i')


subplot(3,2,2)
plot(ss,xx,'LineWidth',1.5)
hold on;
xlabel('s')
title('x')








subplot(3,2,3)
plot(xx,qB,'r--')
 hold on;

plot(xx,qG,'k:')
 hold on;
plot(xx,qs)
 hold on;


subplot(3,2,4)
plot(ss,lambda,'LineWidth',1.5)
hold on;
xlabel('s')
title('\lambda')



subplot(3,2,5)
plot(ss,bx,'LineWidth',1.5)
hold on;
xlabel('s')
title('b')

subplot(3,2,6)
plot(ss,dbxnu,'LineWidth',1.5)
hold on;
plot(ss,dbx,'r--','LineWidth',1.5)
hold on;
% plot(ss,dbx2,'k:','LineWidth',1.5)
% hold on;
xlabel('s')
title('b\prime(s)')


save('datastaticB')
computedynamicB;



print -depsc attemp1_delta08
% basic_static;